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ABSTRACT 


An implicit delta form finite-difference algorithm for Euler 
equations in conservation law form has been used in preliminary 
calculations of three-dimensional wing-vortex interaction. Both 
steady and unsteady transonic flow wing-vortex interactions are 
computed. The computations themselves are meant to guide upcoming 
wind tunnel experiments of the same flow field. Various modifications 
to the numerical method that are intended to improve computational 
efficiency are also described and tested in both two- and three- 
dimensions. Combination of these methods can reduce the overall 
computational time by a factor of 4. 
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The work reported in this document was performed for NASA-Ames 
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I. INTRODUCTION 


Wing-vortex interaction is encountered in many practical applications 
in aerodynamics, and yet it is not a well understood mechanism. Helicopter 
rotor flow fields are a potentially rich source of vortical flow phenomenon. 
The interaction of a trailing vortex wake with the oncoming blades can 
induce unsteady blade loading and aerodynamic noise. Interest in. improving 
helicopter flight performance has provided impetus to develop computer 
simulations to understand the physics of this problem. 

Current numerical algorithms to compute unsteady transonic inviscid 
vortical flow about complex configurations are frequently either inadequate 
or too costly to use for routine design analysis of a large class of two- 
and three-dimensional flow fields. Unsteady potential theory cannot be 
satisfactorialy used for such analyses unless rotational flow effects are 
sufficiently weak and can be modeled with various circulation sheets. 
Numerical algorithms based on the Euler equations are suitable for any 
inviscid flow field simulation, but current numerical algorithms for 
the Euler equations have large computer time and computer storage 
requirements. 

The present study has a two-fold objective. The first of these is 
to apply an existing computer code for computing rotational compressible 
flow field of simple wing-vortex interactions so as to better understand 
the flow phenomenon. These computations are expected to guide follow-on 
wind tunnel experiments and more efficient approximate methods. The 
second objective is to further the methodology and efficiency of current 



numerical procedures. Then, advanced simulations of helicopter flow 
fields can begin as more powerful computers become available. 

With the above objectives in mind, an implicit finite-difference 
procedure for solving the unsteady conservation law equations of inviscid 
flow was upgraded and applied to a simple three-dimensional wing-vortex 
interaction problem shown in Figure 1. Experimental data will be taken 
at a later date for a rectangular wing spanning the wind tunnel walls. 
Interacting with this wing will be an upstream vortex generated from a 
lifting half span wing. The preliminary computations presented here 
model this configuration, although the upstream vortex is analytically 
specified at this time. 

As noted earlier, computations of three-dimensional flow fields using 
Euler equations generally require too much computer time. As a consequence, 
methods of saving computer time and distributing the grid points efficiently 
were studied simultaneously with the wing-vortex interaction calculations. 
Several methods of improving the computational efficiency of the 
present implicit code are demonstrated for two-dimensional 
flows and some are implemented in the three-dimensional code as well. 

In this report the governing equations are reviewed in Section 2 and the 
boundary conditions are discussed in Section 3. The numerical algorithm 
is described in Section 4, with' the relevant grid system in Section 5. 
Results are discussed in Section 6 with concluding remarks in Section 7. 



2. GOVERNING EQUATIONS 


For a perfect gas neglecting viscous effects, the governing partial 
differential equations are the unsteady Euler equations. These can be 
written in a Cartesian coordinates frame of reference and in strong 
conservation law form for three-dimensions as follows: 

‘’t ^ ^ ‘'y ^ ° (1)' 
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The primitive variables of Equation (1) are the density p, the 
three mass fluxes pu, pv, pw in the three coordinate directions 
X, y, and z, and the total energy per unit volume e. In Equation (2), 
p represents the pressure and is nondimensionlized by 'p^; density p by p^; 



u, V, and w the velocity components in x, y, and z directions by a^^^l 
where a^ is the free stream sound speed and y is the ratio of specific 
heats; and e by p a^. 

The pressure, density, and velocity components are related to the 
energy by the equation of state which is written for a perfect gas as 


P ^ /U^ + v^ + W^\ 

( 2 > 


(3) 


In order to transform the physical flow domain into a cubical 
computational domain, the following independent variable transformation 
is employed: 

T = t 

C = s(x,y,z,t) 

n = n(x,y,z,t) (4) 

t = c(x,y,z,t) 


This transformation maps the body and outer boundary surfaces onto 
constant coordinate planes as shown in Figure 2 and thus facilitiates 
the application of boundary condition procedure. The above transformation 
also permits the clustering of grid points in the vicinity of the body and 
in the regions of steep gradients in the flow. 

Application of coordinate transformation. Equation (4), into 
Equation (1) yields the following governing partial differential equations: 
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where 

q = q/0 

E = (5^q + + ?yF + 52S)/J 

F = (n^q + iij^E + tiyF + ti2G)/J (6) 

G = (?^q + + ?yF + 


The flux vectors E, F, and § all can also be written as 



where U, V, and W are the contravariant velocities (without metric 
normalization) given by 


U = 5t + + y + 

W = ?^ + ?j^u + CyV + c^w 
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The metrics required by Equation (7) are obtained from a chain 


rule expansion of x , y , etc. and solved for i y 5 , etc., to yield 
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where O’^ = x^^z^ + y^z^x^ e z^x^^ - x^^z^ - x^^z^ - x^^z^ 
is the transformation Jacobian. 

The transformed equations given by Equation (5) are not much more 
complex than the original set cast in Cartesian coordinates (Eq. (1)). 
Strong conservation-law form of the equations is maintained for shock- 
capturing purposes. 



3. BOUNDARY CONDITIONS 


The tangency condition along the surface c(x,y,z,t) » constant 
Is that W = 0 and Is used In 
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to obtain u, v, and w; 

The pressure along the body surface can be obtained from a normal 
momentum relation, found by combining the three transformed momentum 
equations 


P 






* V'y * 'z'z '’’c 

+ (ixCx * "y'y * 
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- pV(c U +CV +cw) 

^ '^x n y n 


( 11 ) 


where n Is the normal direction to the body surface. 

At the far field boundaries free- stream values are specified except 
on the upstream and downstream grid boundaries. At the upstream grid 
boundary, an analytical vortex Is Initialized on the wing axis and at 
the downstream boundary a simple linear extrapolation is used for all the 
flow quantities. 
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To simulate wind tunnel test conditions at the ends of the wing 
(wing tips), spanwise velocity v is set to zero there and values of 
other flow quantities are set equal to the corresponding values at a 
previous spanwise station. 
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4. NUMERICAL ALGORITHM 


An implicit numerical algorithm is used to solve the equations. This 
method is chosen to avoid restrictive stability conditions when a fine 
grid spacing is used to resolve flow quantities in a region of large 
special flow gradient. * 

4.1. Conventional Form 

The basic numerical algorithm used was developed by Beam and 
1-3 

Warming and has been used extensively and successfully for a large 
class of problems.^ It is second-order accurate in space and time, 
is noniterative, and is in a spatially factored form referred to as the 
"delta- form.'' For either trapezoidal or Euler temporal implicit differenci 
the delta form scheme is given by® 

(I + h6.A" - eTJ""7^Aj)(I + h6 A J)x 

5 ^55 n Inn 

(I + hJ^C" - - 9") = + + j^g") 

• ( 12 ) 

V. ’ 

where the 6's are central difference operators, A and 7 are forward or 
backward difference operators, e.g., A^q = q(5, n, c + A?) - q(§, n, 5). 
Indices denoting spatial location are suppressed for convenience and 
h = At corresponds to Euler implicit first-order and h = At/2 to 
trapezoidal second-order time accuracy. Also, I is the identity matrix 
and^*^* q(nAt) where n is the number of time steps. 



The Jacobian matrices A’’, and are obtained in the time linear- 
ization of and g” and can be written as 
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where 


<1.^ = 0.5(y-1)(u 2 + + w2), 

e = KjU + KjV + K3W 

and for example, to obtain A, 
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Fourth-order dissipation terms such as ^ in 

Equation (12) are added explicitly and these help to control possible 
numerical instabilities. The addition of the implicit second difference 
terms, with coefficient ej, operating on - q*') extends the linear 
stability bound of the fourth-order terms. For more details see 
Reference 6. 

The solution algorithm for the approximately factored form of the 
implicit scheme of Equation (12) takes the following form: 

-(I + hS^A" - ejJ"\^A^J)q = RHS | 

(I + h6^§" - ejJ"V^A^J)q = q \ (14) 

(I + hfi^e" - ejJ"^^A^J)(q'’‘^^ . ^") » q j 

and finally 

-n+1 _ qn ^ ^-n updating the solution 

The central difference operators on the implicit side of Equation (12) 
or (14) produce block-tridiagonal matrix operators which must be inverted 
sequentially to obtain Aq” = (q”"**^ - q”). 

In implementing boundary conditions, the unknown values of q on 
the boundaries are updated explicitly, and Aq is set to zero leading 
to first-order error in time at the boundaries. Explicit treatment of 
the boundaries leads to a simple and flexible scheme, where boundary 
conditions become a modular element which can be put- in or pulled-out 
of a computer code without disturbing the implicit algorithm. 
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The values of p, U, and V along the body surface are found by linear 
extrapolation from above and values of u, v, and w are obtained from 
Equation (10). The surface pressure is obtained by integrating Equation (11). 
In Equation (11), the right-hand side is known from the above extrapolation 
process, and the basic approximate factorization algorithm is applied 
along the body using backward differencing in z and central differencing 
in 5 and n. Scalar tri diagonals are thus inverted in the c and n directions. 

4.2. Diagonal Form 

The solution process of implicit algorithm of Equation (12) involves 
the block-tri diagonal matrix solutions and this constitutes the major part 
of the numerical work. Equations (12) are a coupled set and thereby produce 
a 5 X 5 block structure for the implicit operators of Equation (12). If 
the operators are factored into five scalar operators, the resulting 
system would be more efficiently solved. This idea is used in the diagonal 
form of the implicit algorithm developed in Reference 12. Here, we shall 
briefly outline the development of this algorithm. More details can be 
found in Reference 12. 

The Jacobian matrices A, B, and C (Equation (13)) have a set of 
eigenvalues and a complete distinct set of eigenvectors. Similarity 
transformations (see Warming, Beam, and Hyett, Reference 14) can be 
used to diagonalize A, B, and C where 

A = . § = ri r' . and C . (15a) 
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0 0 U 
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= D[V,V,V,V + c(nj + n^ + n^)*, V - c(nj + nj + np^ ] 
A^ * D[W,W,W,W + c(cJ + cJ + Cp^ , W-c(cJ + cJ + 5^* ] 


(15b) 


where 9 = k u + k v + k w and, for example 
A y i 




(16) 


with k = c for A, k = n for B, and k = ? for C, and 
a = p/(/2c), 6 = l/(/2pc). 

Usinn the similarity transformations. Equation (15a), for the 
Jacobian matrices in Equation (12), the diagonal form of the finite- 
difference algorithm is written in approximate factored form as 
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The new implicit operators (I + etc., of the 

Equation (17) are still bl ock-tridi agonal , but now the blocks are simply 
in diagonal form so that the operators reduce to five independent scalar 
tridiagonal operators. This has a large positive effect on saving the 
overall computational time. 

The solution process for the implicit part of Equation (17) then 
consists of the following steps: 

1. Sj = (T“^)” a matrix-vector multiplication at each grid 
point, since T^^ is known analytically; 

2. Sj = [I + five scalar tridiagonal 

inversions for the operators; 

3. a matrix-vector multiply at each, point; 

4. T = [I + h6 A - a J]“^ S,t five more scalar 

tridiagonal inversions; 

5. = P"^ "S^, a matrix-vector multiply at each grid point; 

6. Sg = [I + h6^A^ - ejO'^v^a^J]”^ five more scalar 

tri diagonal inversions; 

7. Aq’’ = t” Sg, another set of matrix-vector multiply; and 
finally, 

8. q”^^ = q" + Aq” to update the solution. 

This solution procedure constracts with the three bl ock-tridi agonal 
inversions required in Equation (12). 

In both conventional and diagonal schemes discussed above, fourth- 
order central difference operators are used to compute convective 
derivatives. Metrics are computed using second-order central difference 
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formulas at interior points and three-point one sided formulas are 
used at the boundaries. For differencing metrics in three-dimensions 
a free-stream error can be introduced if this differencing is not done 
with special weighted averages.® Perfect maintenance of the free-stream 
can be achieved by simply subtracting the free-stream fluxes from the 
governing equations; that is 

3 q + 3p(^ - E ) + 3 (F - F ) + 3 (S - 6 ) = 0 . 

The metrics are not stored in the present computations because of computer 
storage limitations and hence have to be computed at each time step. 

Equations (12) and (17) can be used for either time-accurate (unsteady) 
or steady. state computations. For steady state calculations, approaches 
zero asymptotically with the solution satisfying the explicit side of the 
equation, which is the exact steady state difference equation. 



5. GRID GENERATION 


Although the transformed equations. Equations (5), are somewhat more 
complicated than the original Cartesian form. Equations (1), they offer 
the following significant advantages: 

(1) The boundary surfaces In the physical plane can be mapped 
on to rectangular surfaces In the transformed plane. 

(2) The grid points can be clustered in regions that experience 
rapid change in the flow field gradients. This Is particularly Important 
in the present problem because of the presence of the Interacting vortex 
and shock waves. 

A surface conforming grid structure Is used to simplify application 
of surface boundary conditions and- to Improve the overall accuracy of 
the numerical scheme. The grid Is generated In two phases. Firstly, an 
0-type, two-dimensional grid Is generated. To take advantage of the 
generality of the transformed equations, one needs a fairly automatic 
method of generating a smoothly varying grid that conforms to arbitrary 
bodies and allows grid point clustering. The elliptic grid generation 
procedure popularized by Thompson, Thames, and Mastin^® with modifications 
by Steger and Sorenson^® and Sorenson^^ was used. Source terms^^ are 
Implemented In the elliptic grid generation equation to cluster points 
both at the body surface and at the center of the grid so as to resolve 
the approaching vortex. Additional source terms^®’^^ are added in the 
present case to pull the grid lines towards the wing axis. The method- 
also enforces orthogonality of mesh lines at the boundaries. Figure 3 
shows several views of a two-dimensional 0-type grid generated by the 
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above method. This grid is then distributed along the rectangular wing 
span using spanwise exponential clustering functions to generate a three- 
dimensional mesh. The grid so generated has 48 points in the periodic 
direction, 20 points in the radial direction, and 21 points along the 
span of the wing. The grid structure is coarse and its dimensions are 
matched to the internal computer storage capability of Control Data 
Corporation 76.00 machine. 



6. RESULTS 


Preliminary numerical experiments using the conventional algorithm 
with the above 48 x 20 x 21 grid for a simple flow past three-dimensional 
wing at = 0.72 and a = 0,deg. revealed that a single solution took 

as much as 8 hours of computer time on CDC 7600 computer. By using the 
diagonal form of the algorithm, this time was reduced by a factor of 
2.5 per time step without the loss of accuracy or numerical stability. 

The magnitude of savings in computational time reported here is much 
better than the previously reported studies with the diagonal algorithm.^^*^^ 
This improvement stems from the fact that in the present application the 
size of the block periodic tridiagonal operator is much more significant. 

The diagonal algorithm is twice as efficient for periodic block tridiagonal 
matrices as it is for a regular block tridiagonal. Because of this big 
savings in computational time, all further computations were undertaken 
using the diagonal algorithm. 

6.1. Two-Dimensional Calculations 

In initial tests of the computer code several two-dimensional flow 
cases were computed. Instead of generating a solution at one spanwise 
station (for a 2-D case), a five plane solution (like a 3-D case) was 
computed with the spanwise index remaining active in the code. The two 
spanwise end planes are updated with simple symmetry conditions. The 
results of such a calculation are shown in Figure 4 in which a plot of 
computed Cp distribution for a steady nonlifting case using the present 
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This solution uses 


code is compared with a reference solution.^® 
the two-dimensional grid shown in Figure 3. Figure 5 shows a C plot, 
for a lifting case. 

In initial calculations, two-dimensional solutions of inviscid free 
shear layers (jets and wakes) with their velocity field superimposed on 
the main stream at the grid boundary of the wing were computed. The excess 
velocity for the jet and the defect for the wake is assumed to be 0.2U 

o» 

on the wing axis. Figures 6 and 7 show the computed steady solution in 
the form of Cp plots. While the shear flows of Figures 6 and 7 are 
symmetric with respect to wing axis. Figure 8 shows a Cp plot for the 
asymmetric case. Asymmetry in the flow is created by locating the axis 
of the symmetric jet of Figure 8 O.OIC above the wing axis. 

6.2. Three-Dimensional Calculations 

As noted earlier, the objective of the present study is to provide 
computational guidance for an upcoming wing-vortex interaction wind 
tunnel experiment. The aerodynamic performance of rotory wing aircraft 
is strongly influenced by the vortex wake shed by the rotating blades. 

The vortex generated by the tip of the blade is time-variant and has 
large strength. The interaction of this vortex with the following rotor 
blade will affect the flow characteristics around the blade and thus 
its aerodynamic loads. In the transonic flow regime this effect is 
more pronounced because of the shock-wave motion. Investigation of the 
wing-vortex interaction is thus important not only for improving the 
understanding of the flow around the wing but also for the optimum design 
of the rotor wing. 



An analytical vortex aligned in the center of the wind tunnel and 
located initially at the grid boundary was made to impinge on the rectangular 
wing as shown in Figure 1. The wing was set at zero angle of attack to a 
free-stream of = 0.72. The specified analytical vortex chosen is a 
steady Lamb or spreading vortex whose cylindrical velocity distribution 
given by^^ 
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2irU r 
00 


(1 - e 


-rVa 
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(19) 


20 

resembles that of an experimentally generated vortex. In Equation (19), 

r 4“ 

^ is the strength of vortex; a = (4vt)* is the core radius; v is the 
kinematic viscosity; t is the time; and r is the radial distance measured 
from the core center. The piressure field induced by this vortex is 
determined using the radial momentum equation 


^ 

dr r 


( 20 ) 


Using the energy equation for constant enthalpy flow. Equation (20) can 
be written as 



( 21 ) 


where is the total enthalpy and q^ = u^ + v^ + w^. The pressure field 
from this and the velocity distribution from Equation (19) are superimposed 
on the main stream at the upstream grid boundary of the wing for initializing 
the vortex. Figure 9 shows a plot of typical radial velocity and pressure 
distribution for a Lamb Vortex of strength = 0.03 at the grid 
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boundary. The three-dimensional computation of wing-vortex interaction 
is carried out using grid dimensions of 48 x 20 x 21 with 21 spanwise 
points. Figure 10 shows the computed steady pressure distribution on the 
wing plotted in the form of Cp distribution at various spanwise stations. 

Because of the coarse grid, it is necessary to consider the vortex 
recovered slightly upstream of the wing leading edge as the input vortex 
for interaction purpose. Typically, the size of this vortex is of the 
order of the span of the wing although the vortex initialized at the grid 
boundary is much smaller. Figure 11 shows a spreading vortex at the grid 
boundary whose radial velocity distribution is given by 


ve(r) 
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2irrU 

00 


(1 - e 


-r^/a^ 



( 22 ) 


Steady state solution of the interaction of this vortex with the wing 
of Figure 1 is computed and a plot of distribution at various spanwise 
stations is shown in Figure 12. Figure 13 shows a detailed Cp distribution 
at a spanwise location of 0.42 chords to the right of the centerline. These 
plots show that the effect of vortex interaction is to superimpose a swirl 
velocity on the free-stream which imparts lift and slight supercritical 
effects . 

Preliminary unsteady computations of interaction is done by choosing 
a vortex with time-dependent strength whose radial velocity distribution 
is given by 
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00 00 

where u is the frequency of the "breathing" motion of the vortex. 

Figure 14 shows a plot of this velocity distribution as a function of uit. 



Interaction of this vortex with the rectangular wing of Figure 1 
is computed for two periods of "breathing" motion at a dimensionless 
frequency (o) of 0.4 and a vortex strength of 0.2. This time period of 
two cycles is considered long enough to depict the unsteady nature of 
the flow field. Figure 15 shows the computed Cp distribution on the 
wing surface at various instants of breathing cycle for this flow field. 

6.3. Further Experiments to Reduce Computer Time 

As indicated previously, computational efficiency is being improved 
by carefully distributing grid points over the body and by using the 
di agonal izati on scheme. Two additional techniques have been demonstrated 
to save overall computational time further and are discussed below. 

(a) Special Variation of Numerical Dissipation Coefficients 

The numerical algorithm, Equation (17), uses added numerical 

dissipation terms to maintain stability for nonlinear equations. So 

as not to degrade solution accuracy, fourth-order dissipation terms are 

added to the explicit side of the equations and second-order terms 

working on - q^) are added to the left hand or implicit side of 

the equations. These terms work well, but analysis of new upwind split 
21 22 

schemes, * which are naturally dissipative, suggest that the 
dissipation terms should be weighted by the "absolute value" of the 

gC 

Jacobian matrices A, B, and C where A = |^ etc. The advantage of using 
such weightings is that dissipation is automatically added as needed for 
variation of the flow parameters and, more crucially, the mesh size 
spacing. 


- 24 - 



« 


Consider the 2-0 system in transformed coordinates 

ii + it + iE » 0 

. 3t 35 3n * 

Using the homogeneous property of flux vectors E(q) and F(q), they 
can be split into two parts via similarity transforms as^^ 


( 24 ) 


A A _ 

E = E ' + E 


where 


and 


E =Aq, E »Aq, 


A = = XAX■^ A* = X ^ ^ 1^1 X’^ 

3^ 2 

Here A is a diagonal matrix whose elements are the eigenvalues of A. 
The eigenvalues of A"*" are nonnegative and those of A“ are nonpositive. 
(In subsonic flow the eigenvalues are of mixed sign since |u| < c.) 

With the use of split flux vectors, one-sided spatial difference 
approximations are possible. For example. Equation (24) can be 
written as^^ 


+A„BJ|,| )]Aq;^|^ 


where 


^ ^j-2 

^5 2A5 

fc = '^^3 ^^.i+1 ' ^.i+2 


6-E = 


2A5 


(25) 
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are backward and forward second-order accurate one-sided difference 
operators, respectively. 

Further, A“ = (A ± |A|)/2 
|A| = X|A|X-" 

= (E ± lA|q)/2 , etc. (27) 

Substituting the relations (27) in (25), we can write Equation (25) as^^ 

[I+At(6^A" + - (2A5 )’^(va)|A1 - (2Ati)"^(va) 1B|] Aq" 

= -At(6^E" + 6 ^f" + (4A5)’'(VA)2|Ajq" + (4An)“"(VA)2|B|q") (2B) 

V + A _ + 5^ 

where 6^ = ^ ^ ^ , etc. 

Calculating the absolute values of the coefficients matrices 
|A| and |B| is very costly. However, the mesh variation enters into 
A and B at a point as a simple scaling. This important variation can 
be accounted for, along with the main flow parameter, by using not [A] 
but its spectral radius (i.e. largest eigenvalue in absolute value). 

Replace |A| with and |Bl with pj^ where p^ and pj^ are the 
spectral radii. For A, the maximum eigenvalue p = |U| + c(c^ + ^5)^* 
and similarly for B, pj^ = 1V| + c(n^ + n^)^. We can further approximate 
these, for computational simplicity, as 

and (29) 

Note that in the standard algorithm, the spectral radius is set to unity. 
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With the above approximations, the new smoothing factor, say in 
the 5-direction, can be written as 

(1+M.)a„ 

'old 5 (l«xl * l«yl> (30) 

where is the old smoothing coefficient. Scaling the dissipation 
terms by the appropriate spectral radius is relatively inexpensive (the 
order of 5X increase in computational work). 

For a three-dimensional test case of a rectangular wing with 
NACA0012 airfoil at o = 0* and = 0.72, a steady solution previously 
obtained in 750 iterations is now obtained'with this scaling in 400 
iterations. This is a savings of over 40% in computer time. The new 
dissipation terms allow use of a larger time step without loss of accuracy, 
as seen in the Cp plot of Figure 16, and thus account for the improved 
efficiency. 

(b) Special Variation of the Time Step 

Keeping in mind that the grid variation accounts for significant 
metric variation, various researchers have long suggested that a 
variable at be used at each point in the computational field. The 
rationale is given that a large At should be used as the mesh cell 
increases in size, or alternatively, the equations should be scaled 
to reduce the condition number. In either event. At should vary as 
the spectral radius of (A6^ + B5^ + C6^), which is an unknown magnitude. 

In our two-dimensional experiments, we have scaled At based on 
simple estimates of the spectral radius of (A6^ + B6^), but we have 



found that a geometric scaling inversly proportional to the Jacobian 
is adequate and far easier to program. The precise algorithm is now 

[(I + 0)h6^A - a»ejJ"^{v^A^)J)(I + 03h6^B - aiejJ‘^(v^A^)J](q"‘^^ - q" ) 

= -o)h(6^&" + - ajE^J’^[(V^A^)^ + (v^A^)^]Jq" (31) 

where u = yj and J = (Sx’^y " results shown in Figures 17 

to 19 indicate the effectiveness of this scaling. For a circular cylinder 
at = 0.45, wave drag can be computed to 4 significant figure accuracy 
within 100 iterations on a 50 x 20 grid using a variable scale factor u 
for the time step At. Using the best constant At (i.e. u>=l), slightly 
over 250 iterations are needed as shown in Figures 17 and 18. 

The circular cylinder calculation uses a stretched polar grid. To 
access the effect of grid distortion, similar runs were made for the 
NACA0012 airfoil at a = 0" and * 0.72 and the grid shown in Figure 3. 
Residual decay using constant At (uj = l) and variation of At is indicated 
in Figure 19 and Figure 20 shows the corresponding Cp distribution. More 
iterations are required for this case (even though it is shock free), than 
the circular cylinder case, but the variable At proves to be at least 
twice as efficient. 

The scaling given above in Equation (31) works effectively when the 
far field grid structure is very coarse, but a scaling such as to = 0”^ is 
unstable as u becomes too large. Figure 21 shows the Cp distribution on 
a circular cylinder at = 0.'45 with a 78 x 26 grid structure computed 
using variable At scheme. 
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7. CONCLUDING REMARKS 


The present study addressed the simulation of transonic three- 
dimensional. wing-vortex interaction. The objective was to guide the 
design of an upcoming wind tunnel experiment and to support the inter- 
pretation and analysis of the data obtained. 

A calculation procedure was adopted which employs an implicit finite- 
difference scheme for solving the Euler equations describing the flow 
about arbitrary geometric shapes. To check the accuracy and resolution 
of the procedure, calculations were carried out for three-dimensional 
wing-vortex interaction, for two-dimensional solutions of simple flow 
past lifting and non-lifting wing as well as for simple shear flows 
superimposed on the main stream. The key to the accuracy of the 
computations was found to be the grid structure. 

The calculations were carried out within the memory of the Control 
Data Corporation 7600 computer. This machine tends to restrict the size 
of the grid which can be used. Consequently, a great deal of effort had 
to be exerted in clustering grid points to flow field gradient regions, 
especially near the vortex core. An elliptic grid generation routine 
with boundary and interior clustering functions was used in this effort 
and greatly aided in the generation of the finite-difference grid. 
Nevertheless, vortices with sharply contained swirl velocity could not 
be adequately resolved with the grid points available. The adequate 
resolution of such regions will likely require the development of more 
sophisticated software such as the overlay of one grid onto another. 
Alternatively, a more advanced computer than the CDC 7600 could be 
employed. 
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Initially computer time also proved to be quite excessive. Typically 
about 10 hours of CDC 7600 time was required to obtain one three-dimensional 
solution. By using a variety of techniques, the overall computational 
time was reduced by a factor of 5. As described in the main text, computer 
time was reduced by using diagonal ization of the implicit operators and 
smoothing terms. Other methods including use of a variable time step 
were tested in two-dimensions but these have not been successfully 
implemented in three-dimensions. It is expected that with further work 
these changes could offer significant improved efficiencies in the current 
and in related computer codes. 

The computations of the wing- vortex. interaction considered a non- 
lifting rectangular wing of NACA0012 airfoil and an analytically represented 
spreading vortex in a free stream of Mach number 0.72. The solutions 
obtained showed that for the vortex considered here, the swirl velocity 
of the vortex imparted lift (rolling moment) and slight supercritical 
effect on the wing. Also, one of the calculations showed that when the 
imposed vortex is large (of the order of the span of the wing), a 
particularly important interaction between the vortex and the tunnel wall 
takes place. This is to be avoided in the design of the experiment. This 
result suggests that if one wants to avoid such an interaction and its 
effect on the flow on the wing, one has to consider a smaller vortex of 
the size of the chord of the wing in the experiment. 

The results reported here are considered as initial in nature and 
are viewed as being qualitatively correct. Any quantitative improvement 
has to come by way of using a much refined grid. 
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Fig. 1. Schematic of experimental arrangement showing 
inertial coordinate system. 



1 . Coordinates transformation and body coordinates. 














= 0.72, a = 0 deg. NACA0012 AIRFOIL 



Fig. 8. Two dimensional asymmetric jet with velocity distribution 

) with 2 = z-0.05 superimposed on 

free stream.. 
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Fig. 9. Velocity and pressure distribution of a Lamb vortex 
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Fig. 14. Velocity distribution of 
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